{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison of `Moran_Local` results with `pygeoda`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from importlib import reload\n",
    "import esda\n",
    "import libpysal as ps\n",
    "import geopandas\n",
    "import pandas\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "PERMS = 9999\n",
    "\n",
    "from libpysal.examples import get_path, load_example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading NCOVR to /home/jovyan/pysal_data/NCOVR\n"
     ]
    }
   ],
   "source": [
    "_ = load_example(\"NCOVR\")\n",
    "\n",
    "nat = geopandas.read_file(get_path(\"NAT.shp\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Load up data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_name = \"HR70\"\n",
    "#y_name = \"HOVAL\"\n",
    "db = geopandas.read_file(get_path(\"NAT.shp\"))\n",
    "\n",
    "w = ps.weights.Queen.from_dataframe(db)\n",
    "cardinalities = pandas.Series(w.cardinalities).values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LISA conditional randomisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `Moran_Local` in `esda`. Single thread."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.01 s ± 11.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "reload(esda.moran)\n",
    "lisa = esda.Moran_Local(db[y_name], w, permutations=PERMS)\n",
    "%timeit lisa = esda.Moran_Local(db[y_name], w, permutations=PERMS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `Moran_Local` in `esda`. Parallel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running parallel implementation on 8 cores\n",
      "249 ms ± 5.66 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "print(f\"Running parallel implementation on {os.cpu_count()} cores\")\n",
    "lisa = esda.Moran_Local(db[y_name], w, permutations=PERMS, \n",
    "                        numba=True, n_jobs=-1)\n",
    "%timeit lisa = esda.Moran_Local(db[y_name], w, permutations=PERMS, \\\n",
    "                                numba=True, n_jobs=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygeoda\n",
    "\n",
    "db_geoda = pygeoda.open(get_path(\"NAT.shp\"))\n",
    "y = db_geoda.GetIntegerCol(y_name)\n",
    "w_geoda = pygeoda.weights.queen(db_geoda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "115 ms ± 925 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "1.1 s ± 3.25 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit lisa_geoda = pygeoda.local_moran(w_geoda, y)\n",
    "lisa_geoda = pygeoda.local_moran(w_geoda, y)\n",
    "lisa_geoda.SetPermutations(PERMS)\n",
    "%timeit lisa_geoda.Run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Comparison of P-values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f07b9a5c490>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfXyU5Z3v8c9vHhIQBSmgCCE8e1JDCdoI4gM+VFt8AhVEVqpnbS3Vs2rr1qdju1rb9bRK1apri9TarqecslZbQMXFrfZUXcUF3UAJqxhRJKRFO8sCozAkk2v/mIRO4pBMMvedebi/79eLVzMzd6657lfqfOe+r+v6XeacQ0REgiuU7w6IiEh+KQhERAJOQSAiEnAKAhGRgFMQiIgEXCTfHeipoUOHujFjxuS7GyIiReX111//s3NuWKbXii4IxowZw7p16/LdDRGRomJmWw/2mm4NiYgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIREQCTkEgIhJwCgIRkYBTEIiIBJyvQWBmM83sLTNrMLNbMrx+mpntMrO6tn+3+dkfEZFCFosneHHzh7y4+QNi8USfva9vJSbMLAw8BJwFNAJrzWylc25Tp0Nfcs6d51c/RETyLRZP0LhzLxWD+wNQ37QbcFSPGMSQQ8uJxRMsfe197v/tZpJtm0aGDRbOGMeVp4xjyKHlvvbPz1pDU4EG59wWADNbBswGOgeBiEjJWlG3nZuf3EDYjH3NSRzQ2vZhHwkZFx03khV120m0dNw2OOngx7/fwk9e2sK986Ywa8pI3/ro562hkcC2tMeNbc91Nt3M1pvZs2ZWnakhM1toZuvMbN2HH37oR19FRDwXiye4+ckN7Gtu5aP9SZLuLyEA0NLqeHxd4ydCIF1LK9z4xAZfbxX5GQSW4bnOZ/sGMNo5VwM8CCzP1JBzbolzrtY5VztsWMYqqiIiBadx517ClumjsGdCZjTu3OtBjw7Svm8tp64ARqU9rgCa0g9wzu12zsXbfl4FRM1sqI99EhHpMxWD+7M/mcy5nZZk8sD4gh/8DIK1wEQzG2tmZcB8YGX6AWY23CwVl2Y2ta0/MR/7JCLim1g8wfpt/3XgNs7LDX/GgxxgwQmjfR0w9m2w2DnXYmbXAKuBMPCoc67ezK5qe30xMBe42sxagL3AfOfcwW+WiYgUqPZB4WgoxP5kki+dPJaH//8WWr1o3OdPRSu2z93a2lqnHcpEpJDE4glO/P7zXQ765iJi8No3z8zpqsDMXnfO1WZ6TSuLRURy9MhLW3wLAYBIOOTrYHHR7VksIlIoYvEEDzz/Nv/46kG3A/aEw/k6WKwgEBHphRV12/nbf6o7sBLYL5EQLJpbU5yDxSIipSoWT3Djr/wLgfKIcc/FUxjYP0r1iIFFXWJCRKQkpNcKGnJoOUtfe5/9HkwLzaQsEmLR3MmcVzPCnzfIQEEgItKF9Gmhza2t/N15x/DQ7xp8ea+ysLHq2pOZcORhvrR/MAoCEZGDSK8VtK9tRcAdT20iGjK8qPwTDRvNSUd52LCQcfecyX0eAqAgEBE5qMade4mGQgdCACBsxt5mb+4L3TdvClXDD+Oj/ckDt53yQUEgInIQFYP709zacW2wVyFw+fTKPh0H6IqCQEQCr/NgcLshh5Zz95zJfONX62nOcYrQIdEQd8+tYV9zkimjDs/LLaCDURCISKB1Hgy+e87kDpvAnDRhKNaLYj9hIP3aoRWYPn5I3m7/dEVBICKBlWkw+IYnNnDMUQMBeLnhQz7ck+hxBdGwwXdmT+K7z2zqEDCFGAKgIBCRAMs0GLy/pZUz73ux120acN8lqa0lZ04anvGWU6FREIhIYFUM7s++Fu9Whp1ZNYy70spBDDm0vKADoJ2CQEQCKRZP8Oo7sZwHgSF1FfCrr55A7dghuXcsDxQEIhIY7bOD1myJcffqt0i2elMsqH80RDQS9qStfFAQiEggLF2zlW8/VY9zjhZPtg37i1bwtUy03xQEIlLylvz+Hf7Ps2961p4BoRCUh0O0QkHPCMqGgkBESk4snqC+aTfgePOPezwLgQFlYZLOcfecyZw0YWhRzAjKhoJARErKirrtfOPxOs9v//SPhrljVjWnVx3RYVZQKdCexSJSMmLxBDc9sd7zEABoda5DCJQSBYGIlIz6pt20+LCJfDRsLJpb3OMAXdGtIREpWAcrBtf5mPqmXbzyToyfvvQuXi0Pu/WcKqqGDwQc1SMGlWwIgIJARApUd8XgAO577i0efKEBL+4EhQALwWUnjObaMyaW9Ad/ZwoCESk4mYrB3fTkBk6aMPTAB/RFP3qZN97f5dl7RiKhvGwTWQgUBCJScDIVg4uGQjTu3AvAd56q9zQEUu0bH/m1I32BUxCISMHJuDNYS5Jn1m/np69s9aw0RLrmZGtRrw7OhWYNiUjBGXJoOfNqKzo815J0LHn5PV9CAOD286sDNS6QTkEgIgWnYccefvlv2/rs/W49u4oFJ4zus/crNL4GgZnNNLO3zKzBzG7p4rjjzSxpZnP97I+IFL6la7Zy1n0velIeuiv9wkY0DHdeMImFp4739b0KnW9jBGYWBh4CzgIagbVmttI5tynDcXcBq/3qi4gUtvTy0N/zsDhcJtGw8e1Z1UwaMagk6gR5wc/B4qlAg3NuC4CZLQNmA5s6HXct8CRwvI99EZECkr5Q7OWGP3PzkxsIAR83+1Abos3Vp45j+vghJb84rDf8DIKRQPpNvkZgWvoBZjYSuBA4gy6CwMwWAgsBKisrPe+oiPSd9IVi+5OttCRb8fMuUDhkrP7aKYFcH5AtP8cILMNznf/cPwRuds51OXnXObfEOVfrnKsdNmyYZx0Ukb6VvlBsT6KFRIv3IWBAeSTEIWVhyiMh7ptXoxDohp9XBI3AqLTHFUBTp2NqgWVmBjAUOMfMWpxzy33sl4jkSePOvURCmb4jeufQ8ggPLTiOQf2jGgPIkp9BsBaYaGZjge3AfODS9AOcc2PbfzaznwNPKwRESkOmgnEbt+8invB39W5zayvVIwYqAHrAtyBwzrWY2TWkZgOFgUedc/VmdlXb64v9em8Rya9MBeNOmjCU21ds9Py9omEjZFAWDh94L4VAz5hz/s7V9Vptba1bt25dvrshIhm0bxH5lcfWkUjbHaZfNMSgfmF27Gn29P3KIyEWzS2tbSP9YmavO+dqM72mWkMi4on2q4AQ1iEEgFQVUY+nhvaLhlhyWS0zjk5NIFEA9J5KTIhIztJnA33c7O0YQHkkxK1nV2V8rXrEQE/fK6h0RSAiOatv2k3Sh8UAYYNFc1Mb0gzoF+GOpzYRDRvJVqexAA8pCESk12LxBEtfe58Hn9+MH4uC759/LOfVjABgwbTRzKwerrEAHygIRKRXVtRt56YnNnxiPMAr0bAxffyQDs8NObRcAeADjRGISI+1jwn4FQKRkHHPxTX60O8juiIQkR5r3LkXr2eeR0KGGVx5yliuPHmcQqAPKQhEJCvpK4V/+vIWz64Gbvz80Uwb+ymikbDu/eeJgkBEutW5YqhXIVAWNk6eOIyaUYd70p70joJARLqUvkZgH96PCQR1w/hCosFiEelSfdNuQhmryucuyBvGFxJdEYjIJzTs2EPdtv8i9tF+7v2Xzb7MDgr6hvGFREEgIh1c9/9eZ+WGP/n6HgPKwkwbN6T7A6VPKAhEBEiNBfzN0tdZ8+5O398r6ZzGBgqIgkBEWFG3neuX1fkwFNzRgLIwSac6QYVGQSAScLF4gpueWO9/CJSHueP8ak6vOkIhUGA0a0gkgGLxBC9u/pAXN3/Aq+/E2N/i3TLhcAi+d+EkyiMdP16SrU4hUKB0RSASEO0rgzdu38XtKzfiR5mgsrDxg4trUmWjyyPc1Gm7SoVAYVIQiJSgzhvHL12zlTue3kTELKeNY4xUVdD9GfYeOHfScL5zwaQDH/azpozUFpJFQkEgUmI6bxw/q+YoHl+3HYD9ObZ91anj+PRRAw9800+0JLn4s6O44qQxTDjysE8cr7LRxUFBIFJCMpWDaA+BXEXDxpWnpKqC6pt+aVEQiJSIWDzB7978gEjI+3IQ0XDH/QH0Tb+0KAhESkD77aAQ8LHHe0aGDZ697pSMt36kNCgIRIpc+u0gr4UN7rtkikKgxCkIRIpcfdNuEj6EQP9IiIcvr2XG0cM8b1sKi4JApEjF4gmWvvY+Dz7/Nh7vGgmAM6geMdCHlqXQKAhEitCKuu3c9MR6Eh6uCAYoC0N5JKIFYAGjIBApMu1jAl6HwGHlER5acCyD+pdpWmjAKAhEikzjzr2+TBFtbm2lesQgBUAA+Vp0zsxmmtlbZtZgZrdkeH22mW0wszozW2dmJ/vZH5FilSoS9wFPr9/OpqZd7MuhTES7M6uOoF80xGHlEfpFQ7oVFGC+XRGYWRh4CDgLaATWmtlK59ymtMOeB1Y655yZTQYeB6r86pNIMVq6Ziu3rdxI0sOJQWGDu+ZOBtAKYfH11tBUoME5twXAzJYBs4EDQeCci6cdPwB8mfwgUrQe/v07fO/ZNz1tM2Jw7yVTOqwSlmDzMwhGAtvSHjcC0zofZGYXAt8DjgDOzdSQmS0EFgJUVlZ63lGRfOhcIbTzaz95aQuLf7/F0/ecUjGIn/718frwlw78DIJMo1mf+MbvnPsN8BszmwF8FzgzwzFLgCUAtbW1umqQote5QujdcyYza8pIIHUr6I6nNrE/x3tBYYPLpo+m4vD+7NiT4AvHHEntWG0YL5/kZxA0AqPSHlcATQc72Dn3opmNN7Ohzrk/+9gvkbzKVCH0pic3cNKEoTz2ynvc/0JDzu8RNlj99RkqDSFZ8TMI1gITzWwssB2YD1yafoCZTQDeaRssPg4oA2I+9kkk7xp37iUaCh0IAYBwyLjqF6+z9r2dObdfFoYfXKz6QJI934LAOddiZtcAq4Ew8Khzrt7Mrmp7fTEwB7jczJqBvcAlzjnd+pGSVjG4P/taOk7//CiR9CQEDikLs/iLxzHj6CNybkuCI6sgMLMTgAeBT5P61h4GPnLOdVmIxDm3CljV6bnFaT/fBdzVwz6LFD2/vu+0Okf1iEG+tC2lK9sFZf8A/BXwNtAfuJJUMIhIDzXu3EtZ2Pu1nGURLQqT3sn61pBzrsHMws65JPAzM3vFx36JlJz26aKvbYl5vnlMWdhYde3JGheQXsk2CD42szKgzszuBv5IagGYiGTByx3EwgYX11bw6ze2Ew2HSDrH3XMmKwSk17INgstI3Ua6Brie1LTQOX51SqRUNOzYw+r6P3Hvc5vJvTpQSjQS4sYvVHHjF6pUHkI8kVUQOOe2tl0RjAF+DbzlnNvvZ8dEit1ty//AY2ve97TNsrB1GAdQAIgXsp01dC6wGHiH1IrhsWb2Vefcs352TqTYtI8DvPvhHu9DIBLSOID4IttbQ/cApzvnGgDMbDzwDKAgEGmTKg1RT0vS4dVQcDRs9IuED5ShUAiIH7INgg/aQ6DNFuADH/ojUnRi8QSPvLSFH3tYIO7rn5vAZdPHACoTLf7rMgjM7KK2H+vNbBWp/QIccDGpEhIigZbaO3gDiRZvp4NWDD5E4wDSZ7q7Ijg/7ecdwKltP38IDPalRyJFIhZP+BICAFNGHe55myIH02UQOOeu6KuOiBSTWDzBnc9s8iUELp9eqbEA6VPZzhoaC1xLavrogd9xzs3yp1sihadhxx7qtv0X2/7zYx78XQOtHpULKguntu64+LOjuOKkMQoB6XPZDhYvB34KPAWeTYgQKRp+rQm44fP/g2njhmgwWPIq2yDY55x7wNeeiBSode/GPA2B6z83gdOqjtSHvxSMbIPgfjO7HXgOSLQ/6Zx7w5deieRR+6KwAWVhVm38Ew88/7ZnbUdC8MXpYxQAUlCyDYLPkKo3dAZ/uTXk2h6LlIz24nCu1ZFIertnQPvOYQoBKTTZBsGFwDjVF5JSlr6XsNe0c5gUsmx3x1gPaGKzlLTGnXtp9mE6KGjnMCls2V4RHAm8aWZr6ThGoOmjUhJi8QSbmnbh1d2gSMgIGZSn1QnSLSEpVNkGwe2+9kIkT2LxBD95aQs/fXkLzV5tGAD0j4Z5aMGxDOpfptlBUvCy3Y/g9353RKSvrajbzvXL6jxZGBOi4wKb5tZWqkcMUgBIUchqjMDM9pjZ7rZ/+8wsaWa7/e6ciB9i8QS/fG0rX/MoBMIG371wEv2iIQ4rj9Avqk3kpbhke0XQYc27mV0ATPWlRyI+aF8bsHH7Lr61fCNeDAWUhQ0zWDS3hllTRjKzerhKRktRynaMoAPn3HIzu8Xrzoj4oX1tQCRkxBPeDQTcO28K08cP6VAuWgEgxSjbonMXpT0MAbXgyZcqEV/5VSo6ajDqU4fog19KQrZXBOn7ErQA7wGzPe+NiEdi8QT1TbtZ8e+NvpSKDoWNisH9PW9XJB+yHSPQvgRSNFbUbecbj9fh1ef//NoKnnij8UB7kVBqXEBXA1Iqutuq8rYuXnbOue963B+RXkkvFHfTE+s9C4Fo2LhxZhU3zqyivmkXYFSPGKgQkJLS3RXBRxmeGwB8GRgCKAgk71L7Bq/HMJqTrTmtDq6tHMT67buJhIxWB4vm/mUaqOoESanqbqvKe9p/NrPDgK8BVwDLgHsO9ntpvzMTuB8IA484577f6fUFwM1tD+PA1c659T05AQm2WDzB3/5TXduHf27zF75QfQQPX3b8gasLTQOVoOh2jMDMPgX8LbAA+EfgOOfczix+Lww8BJwFNAJrzWylc25T2mHvAqc653aa2dnAEmBaz09DgigWT3DV/33dk/pAMyYO4eHLjgc0DVSCp7sxgkXARaQ+oD/jnIv3oO2pQINzbktbW8tIzTQ6EATOuVfSjl8DVPSgfQmoWDzB0tfe54HfbqbFgxAIh+C+S47NvSGRItXdFcE3SFUb/RbwTTNrf95IDRYP7OJ3RwLb0h430vW3/S8Dz2Z6wcwWAgsBKisru+mylLL28YCEFwnQ5uaZVboCkEDrbowg2/0KMrEMz2X8r9fMTicVBCcfpB9LSF2VUFtbq4VsAdS+LiC1OMy7/wvMq61g4YzxnrUnUox6VWIiS43AqLTHFUBT54PMbDLwCHC2cy7mY3+kSLVfBSRbnSfTQr904miGHFrOF6qHM+HIw7r/BZES52cQrAUmmtlYYDswH7g0/QAzqwR+DVzmnNvsY1+kyLTP3GluSXpWKhrgzgsmseCE0R61JlIafAsC51yLmV0DrCY1ffRR51y9mV3V9vpi4DZS6xF+1Db+0OKcq/WrT1IcDmwg7/CsPETYjO9cUM2CaQoBkc7MueK65V5bW+vWrVuX726ITxp27OHs+1/Ey/3jo2Hj2etO0W0gCTQze/1gX7T9vDUkkrX2KaE//O1mWj36bnJIWYhWB3fPmawQEOmCgkDyqmHHHh791/d4fO37nqwJaHfr2VVMGzdEq4NFsqAgkLy5bfkfeGzN+563e/Wp41h4qqaEimQrl3UCIr22/I1tvoRAeSTElaeM87xdkVKmKwLpc5c8/AqvvdttuaqshNuWLR5SFqG5tVWbxov0goJA+sy6d2Pcueo/+PdtuzxpLxIy/vlrpzB4QJmqhYrkQEEgfeKiH73MG+97EwDt7phdfWA2kAJApPcUBOKL9Jr+Nzy+3rMQKI8Yrc7x7fMnaXGYiEcUBOK5pWu28u2nNmIO9nu0MGxebQU3z6zSLSARHygIxFNL12zlm8s3etbelSePZf7xo3QLSMRHCgLxTCye4PaV9Z61FwnB1aeN14e/iM8UBNJr6eMAAHMXv0KLB/UhQkA0YiyaW6MQEOkDCgLplfYKodFQiHiiJadt48MGIYN5tZVceOwIopGwxgFE+pCCQHosFk9w85Mb2Nfcyr4cdwr4ysljOa9mhD74RfJIQSBZSb8NVN+0i4QHdaLP/PQRfPO8YzzonYjkQkEg3WrfKhIHiaQ3JULLIyHumjPZk7ZEJDcKAulSLJ7gG4/XebJXcLvyiLFormoCiRQKBYF0qb5pt2chEA0b150xkUunVSoERAqIgkC6tHvv/px+/4IpR3Hmp49kYP8o1SMGKQBECpCCQDKKxRMs+uc3+ad1jb1uY1bNcH44/zgPeyUiflAQSAexeIIHn3+bn7+6Nad2/vrE0Xx71iSPeiUiflIQCA079rC6/k/8R9Munt64I+f2yiPGtWdM9KBnItIXFAQB5+W+wf2iqZ1PtUuYSHFREARULJ7g1XdinoTASeM/xR2zJvHR/qRWCIsUIQVBAC1ds5XbVmzEg/pwANwxa9KBMtEiUnwUBAGQXh7isVfe4/4XGjxr+/LplQoBkSKnIChxqfIQGwgZ7PWgPlC7aNj45ZXTqB07xLM2RSQ/FAQlLBZPcMOv1tPsUX0gA/qXhWh1qQFhhYBIaVAQlKiGHXv4+SvveRYCh0RDLL6slkH9oxoQFikxIT8bN7OZZvaWmTWY2S0ZXq8ys1fNLGFmN/jZlyC5bfkfOPO+F/nFa95MCwVoBapHDKRm1OEKAZES49sVgZmFgYeAs4BGYK2ZrXTObUo77D+B64AL/OpHULQPCDe3JD1bFxA26BcNk3ROawNESpift4amAg3OuS0AZrYMmA0cCALn3AfAB2Z2ro/9KHnp20bubUnm3F7VkQNY+pXpAAdmGykEREqXn0EwEtiW9rgRmNabhsxsIbAQoLKyMveelRAvt40EOOPoYTz6pakHHisAREqfn2MEluG5Xo1cOueWOOdqnXO1w4YNy7FbpaW+abcn20aawa3nVHUIAREJBj+vCBqBUWmPK4AmH98vMNrHA55e38RPXn43p7YuOb6Ccz8zguoRA/XtXySg/AyCtcBEMxsLbAfmA5f6+H6B0L5ArDnZmnOJiAtqjuKuOTXedExEipZvQeCcazGza4DVQBh41DlXb2ZXtb2+2MyGA+uAgUCrmX0dOMY5t9uvfhUzrxeI/d351Z60IyLFzdcFZc65VcCqTs8tTvv5T6RuGUk31r0b45GX3805BEJAOAT3zJuiW0EiAmhlccGLxRPMe/gV3vnw45zamTDsEBZ/sValokXkExQEBah9MPiF/9iRc6VQAx65/LN87pjh3nROREqOgqCAxOIJHnj+bX6xZiutrpdzbdNEQnDvvCkKARHpkoKgQCxds5W/82izmM9VDeN/njhWU0JFJCsKggKwdM1Wvrl8oydtDSgLc93njqZm1OGetCcipU9BkEexeIL6pl3ctsKbEABIOkfF4P6etScipU9BkCftheKakw6PlgVQHjFVCRWRHlMQ5EHDjj3c8Ph6mnMcEAgZ9C8L05J0XHP6BC6dVqkQEJEeUxD0sfueeyvnKaE/mDuZKaMOZ/CAMpWJFpGcKQj6SCyeYM6P/5X3Yntzauf1b53Z4UNfASAiuVIQ+CwWT7D0tfe5918259ROVGUhRMQnCgIftVcKTbTktl/A1aeO48pTxikERMQXCgIfxOIJfvNGI3c++yaul+PB/SIhWoHbzz+GBdNGe9o/EZF0CgIPtd8Guu9fNudUHuJznx7GdWccrUFgEekTCgKP3PfcW/zD7xpyXhNQFg5x95waBYCI9BkFQY5i8QSX/uRV3trxUc5tRcPGDy7WgjAR6VsKghwsXbOVby3fmHOV0BDwwF8dy/TxQxQCItLnFAS9EIsnePD5t/n5q1tzaicERCPGork1nFczwpvOiYj0kIKgh1bUbef6ZXXkMiHUgP99dhXTxg3RgLCI5J2CoAee3/QnvrasLqc2hg6IsPr60/ThLyIFQ0GQpfMfeJE/NO3p9e8b8JVTxnDrudXedUpExAMKgm7E4glOW/QCexK9vxn01yeO5tozJuoqQEQKkoLgIGLxBDc+sZ4X3vyw121MHT2YH1/2WQWAiBQ0BUEGK+q25zwWUB4JKQREpCgoCDpZ/sY2vv74hl7//mHlEZpbW7VTmIgUDQVBmppv/zO79iV79bunHz2MH8yr0UYxIlJ0FATAundjzH14Ta9//9Zzqlg4YzygjWJEpPgEPgimfvc5Pviouce/d0gEFl8+leoRA/XhLyJFLbBBEIsn+Ozf/7ZXv3v60cP42ZemetwjEZH8CPnZuJnNNLO3zKzBzG7J8LqZ2QNtr28ws+P87E+7FXXbexUCBvz2+hkKAREpKb5dEZhZGHgIOAtoBNaa2Urn3Ka0w84GJrb9mwb8uO1/fROLJ3o1NfQzIw/jqWtn+NAjEZH88vPW0FSgwTm3BcDMlgGzgfQgmA085pxzwBozO9zMjnLO/dGvTvX0SiAEPP7VE6gdO8SfDomI5JmfQTAS2Jb2uJFPftvPdMxIoEMQmNlCYCFAZWWl5x09mL+fdQxfPHFsn72fiEg++BkEluG5znu4ZHMMzrklwBKA2traXPeB6dZXThrNN8+f5PfbiIgUBD+DoBEYlfa4AmjqxTGeeu/75zLmlme6fF1EJEj8nDW0FphoZmPNrAyYD6zsdMxK4PK22UMnALv8HB9ol+nD/ofzJisERCSQfLsicM61mNk1wGogDDzqnKs3s6vaXl8MrALOARqAj4Er/OpPZ/rQFxFJ8XVBmXNuFakP+/TnFqf97IC/8bMPIiLSNV8XlImISOFTEIiIBJyCQEQk4BQEIiIBpyAQEQk4BYGISMApCEREAk5BICIScJZa01U8zOxDYKtHzQ0F/uxRW8UgSOercy1dQTpfL891tHNuWKYXii4IvGRm65xztfnuR18J0vnqXEtXkM63r85Vt4ZERAJOQSAiEnBBD4Il+e5AHwvS+epcS1eQzrdPzjXQYwQiIqIrAhGRwFMQiIgEXCCCwMxmmtlbZtZgZrdkeN3M7IG21zeY2XH56KcXsjjXKjN71cwSZnZDPvropSzOd0Hb33SDmb1iZjX56KcXsjjX2W3nWWdm68zs5Hz00wvdnWvaccebWdLM5vZl/7yWxd/2NDPb1fa3rTOz2zztgHOupP+R2ibzHWAcUAasB47pdMw5wLOAAScAr+W73z6e6xHA8cCdwA357nMfnO+JwOC2n88u8b/tofxl3G8y8Ga+++3XuaYd9wKpXRDn5rvfPv9tTwOe9qsPQbgimAo0OOe2OOf2A60DiQ8AAANySURBVMuA2Z2OmQ085lLWAIeb2VF93VEPdHuuzrkPnHNrgeZ8dNBj2ZzvK865nW0P1wAVfdxHr2RzrnHX9qkBDACKdSZINv/NAlwLPAl80Jed80G25+ubIATBSGBb2uPGtud6ekwxKJXzyFZPz/fLpK78ilFW52pmF5rZm8AzwJf6qG9e6/ZczWwkcCGwmOKX7f+Pp5vZejN71syqvexAEILAMjzX+ZtSNscUg1I5j2xlfb5mdjqpILjZ1x75J6tzdc79xjlXBVwAfNf3Xvkjm3P9IXCzcy7ZB/3xWzbn+wapWkE1wIPAci87EIQgaARGpT2uAJp6cUwxKJXzyFZW52tmk4FHgNnOuVgf9c1rPfrbOudeBMab2VC/O+aDbM61FlhmZu8Bc4EfmdkFfdM9z3V7vs653c65eNvPq4Col3/bIATBWmCimY01szJgPrCy0zErgcvbZg+dAOxyzv2xrzvqgWzOtZR0e75mVgn8GrjMObc5D330SjbnOsHMrO3n40gNPBZj8HV7rs65sc65Mc65McATwP9yznn6LbkPZfO3HZ72t51K6rPbs79txKuGCpVzrsXMrgFWkxqdf9Q5V29mV7W9vpjUrINzgAbgY+CKfPU3F9mcq5kNB9YBA4FWM/s6qRkKu/PW8V7K8m97GzCE1DdGgBZXhJUrszzXOaS+0DQDe4FL0gaPi0aW51oysjzfucDVZtZC6m8738u/rUpMiIgEXBBuDYmISBcUBCIiAacgEBEJOAWBiEjAKQhERAJOQSDSQ2ZWYWYrzOxtM3vHzO43s7K2CpFPH+R33ivSxV0SAAoCkR5oW9Tza2C5c24icDSpqp935rVjIjko+QVlIh47A9jnnPsZgHMuaWbXA+8Cv2s/yMyGAL8EhgH/RuZ6MiIFQVcEIj1TDbye/kTbquz3gQlpT98OvOycO5ZUuYDKPuuhSA/pikCkZ4zMFU47Pz8DuAjAOfeMme3M8DsiBUFXBCI9U0+q8uUBZjaQVPXIdzodq/otUhQUBCI98zxwiJldDmBmYeAe4OekCha2exFY0HbM2cDgvu2mSPYUBCI90Fbx8ULgYjN7G9gM7ANu7XToHcAMM3sD+DypMQSRgqTqoyIiAacrAhGRgFMQiIgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIREQC7r8BV6rC0Q3CZB4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pandas.DataFrame({\"Old\": lisa_old.p_sim,\n",
    "                  \"Numba\": lisa.p_sim\n",
    "                 }).plot.scatter(\"Old\", \"Numba\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
